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Abstract 


This  report  is  the  first  in  a  series  presenting  a  study  of  self-correlation 
algorithms  in  intelligence  systems.  It  was  performed  by  the  Mathematics 
Clinic  of  the  Claremont  Graduate  School  in  support  of  the  Algorithm  Analysis 
subtask  of  the  U.S.  Army  Intelligence  Center  and  School  (USAICS)  Software 
Analysis  and  Management  System  (USAMS)  task  at  Jet  Propulsion  Laboratory.  '  The 
self-correlation  algorithms  use  multivariate  statistical  tests  to  determine 
the  equality  of  mean  vectors  from  two  different  datasets.  The  statistical 
tests  developed  were  variations  of  Hotelling's  T^-statistics.  The  main  results 
deal  with  the  analysis  of  the  robustness  of  these  statistics  with  respect  to 
normality  and  equal  covariance  matrices.  To  do  this  the  Clinic  developed 
mathematical  techniques  and  simulation  programs  which  generate  multivariate 
normal,  gamma,  and  lognormal  distributions  with  a  specified  dependency  between 
the  vector  components.  Theoretical  and  sample  multivariate  skewness  numbers 
of  these  distributions  were  also  developed.  The  results  indicate  that  the 
statistical  tests  are  sensitive  to  skewness  in  the  distributions  but  are  not 
affected  by  slightly  differing  covariance  matrices. 


Key  Words: 


Self-correlation  algorithm 
Multivariate  distributions 
Hotelling's  T^-statistic 


i- 


Robustness 

Multivariate  skewness 


„y  This  study  is  the  first  in  a  series  of  reports  involved  in  researching 

self-correlation  and  cross-correlation  algorithms  in  intelligence  systems. 

These  algorithms  are  used  to  maintain  a  data  base  of  current  information  about 

a  battlefield.  The  initial  view  of  the  battlefield  is  stored  in  a  central 

computer  data  base.  As  new  data  is  received  from  the  sensors  on  the 

battlefield,  it  is  used  to  update  the  old  data  and  formulate  a  new  picture  of 

j*£  the  battlefield.  The  work  on  these  algorithms  reported  here  focuses 

on  the  sensitivity  of  the  mathematical  tests  to  changes  and  uncertainties  r(.  , 

-  1 

in  the  data.  J  This  project  was  performed  by  the  Mathematics  Clinic  of  the 

Claremont  Graduate  School  (CGS)  in  support  of  the  Algorithm  Analysis  subtask 

of  the  U.S.  Army  Intelligence  Center  and  School  (USAICS)  Software  Analysis 

tm 

and  Management  System  (USAMS)  task  at  Jet  Propulsion  Laboratory  (JPL).  This 
is  an  ongoing  task  to  study  intelligence  algorithms  for  the  Combat  Developers 
Support  Facility  at  USAICS. 

The  analysis  this  year  focused  on  statistical  tests  that  are  variations 

2  2 
of  Hotelling's  T  -statistic.  In  many  applications  the  Hotelling's  T  -statistics 

2  2 

revert  to  chi-square  (x  )  tests.  (The  x  forms  are  used  in  many  of  the 
current  systems.)  These  tests  are  used  to  check  the  equality  of  means  in  a 
multivariate  setting.  Such  tests  could  be  used,  for  example,  to  test  if  two 
location  means  or  if  two  radar  signal  means  are  the  same.  The  statistical 
tests  have  several  variations  depending  on  whether  the  mean  of  the  old 
data,  and  the  dependency  relationships  of  the  old  data  and  new  data  are 
assumed  to  be  known  or  are  estimated.  Three  main  cases  were  studied: 


1.  The  mean  of  the  old  data  is  assumed  to  be  known.  The 
dependency  relationships  for  both  the  old  and  new  data  are 
assumed  to  be  known. 

2.  The  mean  of  the  old  data  is  assumed  to  be  known.  The  dependency 
relationship  is  estimated  for  the  new  data  and  is  assumed  to  be 
known  for  the  old  data. 

3.  The  mean  of  the  old  data  is  estimated.  The  dependency  relationships 
for  both  the  old  and  new  data  are  estimated  and  are  assumed  to  be 
equal. 

Once  the  statistical  tests  were  developed,  the  Clinic  team  began  investigating 
the  robustness  of  the  tests,  i.e.  the  sensitivity  of  the  tests  to  the 
relaxation  of  the  assumptions.  All  of  the  tests  assume  that  the  incoming  data 
is  normally  distributed  (Figure  i).  This  assumtion  was  deemed  most  likely 
to  fail  in  two  ways:  the  distributions  are  skewed  or  they  have  fat  tails.  To 
test  the  behavior  of  the  tests  when  the  distribution  of  the  incoming  data  was 
skewed,  the  gamma  (Figure  ii)  and  lognormal  (Figure  iii)  distributions  were 
used.  These  distributions  were  chosen  because  they  are  representative 

* 

of  skewed  distributions  and  have  statistical  properties  which  made  the 
mathematical  development  possible.  A  major  task  of  the  Clinic  team  this  year 
was  to  develop  a  software  package  to  simulate  representative  multivariate 
data  from  these  distributions.  The  results  obtained  indicate  that  if  the 
actual  distribution  of  the  data  is  skewed,  the  statistical  tests,  developed 
under  the  normality  assumptions,  may  be  invalid.  The  simulation  programs  were  als 
also  used  to  test  the  robustness  of  the  assumption  of  equal  covariance  matrices 
while  the  investigation  of  fat-tailed  distributions  was  left  to  future  studies. 

The  development  of  simulation  programs  for  determining  the  robustness  of 
the  intelligence  algorithms  and  the  analysis  based  on  the  statistical  tests 
is  being  continued. 
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Chapter  I  -  Introduction 


PREVIOUS  PAGE 
IS  BLANK 


The  Mathematics  Clinic  project  involved  researching  methods  to  maintain 


a  data  base  of  current  information  about  a  battlefield.  The  initial  view  of 
the  battlefield  is  stored  in  a  central  computer  data  base.  As  new  data  is 
received  from  the  sensors  on  the  battlefield,  it  is  used  to  update  the  data 
base  and  formulate  a  new  picture  of  the  battlefield. 

The  process  can  be  illustrated  by  Figures  1  and  2.  The  original  view 
of  the  battlefield,  illustrated  by  Figure  1,  shows  one  plane,  a  truck  with  two 
radars,  and  a  single  radar  on  top  of  the  hill.  New  data  is  collected  from 
the  battlefield,  as  represented  by  the  dotted  lines  in  Figure  2.  It  indicates 
that  the  previously  identified  truck  and  plane  have  moved  to  new  locations, 
and  that  the  truck  now  carries  only  one  radar.  It  also  shows  a  new  incoming 
plane.  The  single  radar  on  top  of  the  hill,  however,  is  still  at  the  same  location. 
The  Clinic's  goal  was  to  develop  statistical  methods  to  analyze  the  incoming 
data  in  order  to  accurately  update  the  data  base. 

The  procedure  followed  to  formulate  the  new  estimate  of  the  battlefield 
is  illustrated  in  Figure  3.  New  data  is  first  collected  and  filtered.  This 
procedure  prepares  the  data  for  the  identification  stage,  where  the  source  of 
the  transmission  can  be  determined.  It  is  the  next  stage,  self-correlation,  on 
which  the  Clinic  team  concentrated  its  research  efforts.  Statistical  tests 
are  performed  here  to  compare  the  new  incoming  data  with  the  old  estimate 
of  the  battlefield.  The  new  estimate  of  the  battlefield  is  then  stored  in  the 
data  base.  Figure  4  shows  the  updated  picture  of  the  battlefield.  The  final 
stage  is  battlefield  identification  or  cross-correlation  which  is  the  analysis 
of  the  configuration  of  the  battlefield  based  on  the  enemies  capabilities. 

The  statistical  tests  used  in  the  self-correlation  stage  are  variations 

of  the  Hotelling's  T  -statistic,  which  is  a  multivariate  extension  of  the 

* 

student  s  t-statistic.  The  use  of  a  multivariate  test  statistic  is  necessary 
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Figure  2 
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Procedure  to  Update  Data 
Figure  3 


because  the  data  from  the  sensors  comes  in  the  form  of  vectors.  (The  location 
vector,  for  example,  may  have  three  components:  longitude,  latitude,  and 
altitude.  The  signal  parameter  vector  may  have  four  components:  pulse 
repetition  interval,  pulse  width,  scan  rate,  and  frequency.)  During  the  first 
semester,  the  Clinic  team  rederived  the  basic  Hotelling's  T2-statistic  and  its 
variations  and  then  extended  the  results.  Chapter  II  presents  the  mathematical 
models  of  the  test  statistics  that  the  team  developed. 

2 

The  statistical  tests  used  in  this  report,  namely  Hotelling's  T  -statistic 
and  its  variations,  were  derived  under  certain  assumptions.  In  some  applications 
these  assumptions  may  not  be  met.  Thus  it  is  desirable  to  investigate 
the  robustness  of  the  statistical  tests  used.  (A  statistical  test  is  robust 
if  the  failure  of  assumptions  does  not  invalidate  its  ability  to  yield  correct 
results).  One  of  the  assumptions  that  was  relaxed  was  that  of  normally 
distributed  data.  For  example,  the  normality  assumptions  are  violated  when  the 
data  comes  from  skewed  distributions  such  as  gamma  or  lognormal.  To  measure 
the  non-normality  of  multivariate  distributions  K.  V.  Mardia's  generalizations 
of  skewness  and  kurtosis  were  used.  Chapter  IV  develops  these  concepts.  A 
second  assumption  which  was  relaxed  was  that  of  equal  covariance  matrices. 

Chapter  V  contains  a  brief  discussion  of  covariance  matrices.  To  study  the 
robustness  of  the  above  two  assumptions,  the  Clinic  team  needed  to  generate 
multivariate  distributions  that  had  characteristics  representative  of  the 
battlefield  data.  This  involved  generating  random  samples  of  vectors  such 
that  the  vectors  components  have  different  statistical  distributions  which 
may  not  be  independent  of  each  other.  That  is  although  the  vectors  in  the 
generated  random  sample  are  independent  of  each  other,  the  statistical 
distributions  of  the  vector  components  may  be  dependent.  The  dependence  among 
the  vector  components  is  described  by  the  covariance  matrix  of  the  multivariate 
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distribution.  For  example,  if  the  covariance  matrix  has  non-zero  off- 
diagonal  elements,  then  a  dependency  between  the  vector  components  exists. 

The  method  devised  for  generating  these  random  samples  consists  of  first 
specifying  the  desired  component  distributions  and  then  specifying  the 
component  dependency  (i.e.  the  covariance  matrix).  The  components  of  the 
sample  vectors  are  then  assumed  to  be  expressable  as  a  combination  of  appropriate! 
chosen  univariate  distributions.  Using  the  assumed  covariance  matrix,  it  is 
sometime  possible  to  solve  for  the  coefficients  of  the  combinations  and 
hence  obtain  the  desired  random  sample.  Thus  the  problem  of  simulating  data 
sets  is  reduced  to  finding  the  coefficient  matrix  of  a  system  of  equations. 

The  Clinic  team  named  the  process  of  finding  the  coefficients  matrix 
"backsolving"  since  the  algorithm  starts  with  the  desired  result  (the 
covariance  matrix)  and  ends  with  the  means  by  which  the  problem  is  solved 
(the  coefficient  matrix).  The  details  of  the  "backsolving"  process  are 
discussed  in  Chapter  III.  The  associated  computer  simulation  programs  are 
described  in  Chapter  VI. 

The  results  of  the  Clinic  team's  investigations  indicate  that  the 
statistical  tests  are  sensitive  to  skewness  in  the  distributions  but  are 
not  affected  by  slightly  differing  covariance  matrices.  A  detailed 
discussion  of  the  robustness  results  are  contained  in  Chapter  VII. 

The  report  concludes  with  Chapter  VIII  which  contains  conments  on 


possible  future  projects. 


Chapter  II  -  Mathematical  Models  for  the  Self-Correlation  Stage 

The  Clinic  Team  analyzed  various  models  for  the  self-correlation  stage. 

In  this  chapter  we  briefly  illustrate  the  self-correlation  algorithm  and 
then  develop  the  appropriate  mathematical  models. 

The  self -correlation  stage  involves  the  testing  for  equality  of  mean 
vectors  from  two  datasets.  The  first  dataset  represents  the  old  estimate  of 
the  battlefield  which  has  been  stored  in  a  central  computer  system  and  is 
referred  to  as  the  Base  Data.  The  second  dataset  represents  the  new  estimate 
of  the  battlefield  and  is  referred  to  as  the  New  Data.  Information  from 
both  datasets  is  summarized  and  stored  in  the  form  of  (mean)  vectors  and 
(covariance)  matrices;  denoted  by  ^  and  z  respectively. (Vectors  will  be 
denoted  by  an  underscore.)  The  covariance  matrix  desribes  the  relationships 
between  the  vector  components.  When  referring  to  the  Base  Data  we  will  use  the 
subscript  B  and  when  referring  to  the  New  Data  we  will  use  the  subscript 
N;  for  example,  ig  denotes  the  covariance  matrix  from  the  Base  Data. 

Figure  5  shows  how  the  self-correlation  algorithm  is  implemented  for  the 
situations  mostly  likely  to  occur  on  a  battlefield.  For  simplicity,  not  every 
possible  case  has  been  included;  however,  the  logic  used  for  most  standard 
situations  is  illustrated.  The  first  step  in  the  algorithm  is  to  test  whether 
the  candidate  entity  (New  Data)  is  in  the  same  location  as  a  known  entity 
(Base  Data).  If  the  locations  of  the  two  entities  are  determined  to  be  the 
same,  then  their  signal  parameters  are  compared.  If  these  two  sets  of 
parameters  are  compatible  with  one  another,  then  we  conclude  that  the 
observations  refer  to  the  same  entity.  However  if  the  candidate  entity 
is  strongly  associated  with  more  than  one  known  entity,  then  it  is  possible 
that  either  the  candidate  entity  or  the  one  of  the  known  entities  is  a 


deceptive  reading  referred  to  as  a  phantom.  Phantoms  can  be  caused  by 
failing  to  divide  the  observations  of  the  two  entities  in  the  separation 
stage  or  by  incorrectly  matching  a  candidate  entity  with  a  known  entity.  When 
a  candidate  entity  has  the  same  location  as  a  known  entity  but  different  and 
noncompatible  signal  parameters,  the  possibility  that  it  is  a  known  entity 
which  has  changed  its  signal  parameters  must  be  checked.  The  last  possibility 
we  consider  is  when  the  candidate  entity's  location  is  not  the  same  as  the 
location  of  any  known  entity.  In  this  case  it  must  be  determined  whether 
the  candidate  entity  is  a  known  entity  with  the  same  signal  parameters  that 
has  moved  (i.e.  the  known  entity  is  mobile)  or  whether  the  candidate  entity 
is  a  previously  undetected  entity. 

The  questions  asked  in  the  above  decision  process,  as  well  as  others  in 

the  self-correlation  algorithm,  are  answered  by  using  hypothesis  tests  .-'f  the 

equality  of  mean  vectors.  The  desired  statistical  tests  are  based  on  Hotelling's 
2 

T  -statistic:  or  a  variation.  In  the  development  of  these  tests  there  are  two 
possibilities  that  must  be  considered;  namely  whether  the  mean  vectors  and 
covariance  matrices  are  to  be  considered  as  known  constants  or  as  estimates. 

In  all  there  are  six  variations  of  the  model  to  be  considered.  These  six 

cases  are  described  verbally  in  Figure  6.  In  Figure  7  the  six  cases  are  enumerated 

and  described  using  previously  introduced  notation. 

As  an  example,  under  the  assumption  that  wB  is  estimated,  is 
unknown  but  Eg  *  E^  is  assumed,  case  5  is  appropriate.  Cases  1-5  lead  to 
variations  of  Hotelling's  T2-statistic  to  test 

V  U*  = 

which  is  called  the  null  hypothesis  against 

HA:  ^ 

which  is  called  the  alternative  hypothesis. 


KNOWN 

BASE  DATA 

ESTIMATED 

f-  ■  - - 

BASE  DATA 

Dependency 
relationship 
of  New  Data 
is  known 

Dependency 
relationships 
for  Base  and  New 
Data  are  known 
and  equal 

Dependency 
relationships  for 

Base  and  New  Data 
are  known  and 
unequal 

Dependency 
relationship 
of  New  Data 
is  not  known 

Dependency 
relationships 
are  not  known, 
but  are  assumed 
to  be  equal 

Dependency 
relationships 
are  not  known  and 
and  are  assumed  to 
to  be  unequal 

Description  of  Mathematical  Models 
Figure  6 
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i>g  known 

Ug  estimated  j 

Case  1 

Case  3 

Case  4 

known 

Z|Y  known 

E^  known 

EB  known 

Eg  known 

eB  =  rN 

E  i  E 

B  *  N 

Case  2 

Case  5 

Case  6 

E^  estimated 

en  estimated 

E^  estimated 

Eg  estimated 

Eg  estimated 

ZB  =  EN 

eB  *  en 

Numeration  of  Mathematical  Models 
Figure  7 
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Before  we  develop  the  mathematical  models  for  these  six  cases,  we  give 


a  more  detailed  description  of  the  mean  vector.  A  sensor  observes  p 
characteristics  of  a  particular  entity  N  times:  therefore,  there  are  N 


vectors  of  the  type  xs 


The  incoming  data  is  assumed  to  be  from  a 


p-variate  normally  distributed  population  and  the  N  observations  are  assumed 
independent  of  one  another.  The  jth  component  of  each  vector  satisfies  the 
equation 

XjU)  -  ,j  ♦  <j(t) 

where  x.(t)  is  the  value  of  the  jth  characteristic  at  time  t,  u-  is  the 
J  J 

true  value  of  the  characteristic,  which  remains  constant  over  time  for  a 
stationary  entity  and  €  . ( t )  is  the  error  term  at  time  t.  The  N  vectors 

•J 

are  averaged  to  obtain  the  sample  mean  vector 


5  (t)  *  £  +  i  (t) 

That  is,  x  (t)  is  the  sample  mean  observed  by  the  sensor  at  time  t  (average 
or  final  time  of  observations).  We  now  proceed  with  development  of  the 
mathematical  models  'or  the  self-correlation  stage 

Case  1  is  the  simplest  and  least  likely  case  to  occur,  the  reader 
Is  referred  to  Johnson  and  Leone  {1964,  pp.  294-295).  The  T2 -statistic 
used  In  Case  2  Is  presented  next. 

Let  Xj , . . . , x^  be  N  observations,  each  with  p  characteristics. 
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The  distribution  of  x^  is  assumed  to  be  normal  with  mean  £  and  covariance 
matrix  £  (In  the  self-correlation  algorithm  £  =  and  z  =  z^.)  Each 

observation  vector  has  probability  density  function 


t  (x^  =  - ! - 

(  /?Tt  )P  |  E  |  ^ 


—  expf-4(x.-  u)'  z_1(x.-  M)] 


The  joint  density,  or  likelihood  function,  is 


V..*  <£!»•••  *2^)  =  L  (£,'  E-1) 


-1  *  ’  ’  ' *— N 


=  11  ^fx 
a=l  ““ 


■1  iN/2 


expL'4  (^,] 


L(£',  £  )  denotes  the  likelihood  function  of  the  parameters  listed.  For 

this  case,  the  null  hypothesis  is  HQ:  u  =  ^  with  ^  known,  and  z  is 
unknown  and  invertible.  (In  the  self-correlation  algorithm,  £  =  yN,z  = 

In  order  to  test  we  compute  the  likelihood  ratio: 


H'“  “N 


max  L(pn,Z_1) 

-1 

z 


max  L(£M_i) 
-1 

v  ,£ 
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The  numerator  of  this  ratio  is  the  maximum  of  the  likelihood  function  with 
parameter  space  restricted  by  the  null  hypothesis,  u».  The  denominator  is  the 
corresponding  function  maximized  over  the  entire  parameter  space,  n.  The 
test  is  to  reject  Hq  if  the  computed  x  is  small;  that  is,  less  than  or 
equal  to  some  Aq.  The  derivation  of  the  maximum  of  the  likelihood  function 
over  the  entire  parameter  space,  fi,  and  over  the  restricted  parameter 
space,  w,  is  shown  in  Appendix  I.  Using  this  derivation  the  likelihood 
ratio  becomes 


|A  +  N  (x  -un)1N/2 


where  A  =  z  (x  -  x)(x  -  x)'  *  (N  -  1)S,  with  S  =  A/(N-1).  Thus 


_  t  —  —a 

a-  1 


.  2/N  _ 


or  equivalently 


I A  +  [/N(x  -  u0)J[/N(x  -  uQ)]- 


1  +  TV(N-l) 


where 


T2  =  (N  -  l)N(x  -  j^'A^lx  -  Uq) 
=  N(x  -  EqJ'S"1 (x  -  u0) 


Equation  2.1  can  also  be  expressed  as 

T2  =  (N-1)(a'2/N-1). 


The  decision  whether  or  not  to  reject  the  null  hypothesis  is  determined 
by  comparing  A  against  a  critical  Aq  at  a  certain  significance  level.  We 
reject  HQ  when 


2/N  2/N 
a  _aq 


(2.2)  - 


Inverting  Equation  2.2  ,  subtracting  1,  and  multiplying  by  N-l  ,  the  critical 
region  is  redefined  as 

2  2 
T  >  T 

i  _  i0 

2  ... 

In  this  case,  T  can  be  shown  to  be  distributed  according  to  an  F  distribution 

with  p  and  N-pH  degrees  of  freedom.  This  result  is  proved  in  Appendix  II- 

2 

For  cases  3  through  6,  the  two  sample  tests,  the  T  -statistic  does  not 
follow  one  specific  distribution;  it  varies  depending  on  the  assumptions 
on  the  covariance  matrices.  When  the  covariance  matrices  are  unknown,  and 

there  is  not  enough  data  to  assume  that  the  estimates  are  equal  to  the  true 
values,  we  use  statistical  methods  to  test  whether  or  not  two  covariance 
matrices  are  equal.  The  derivation  of  this  test  is  presented  in  Appendix  III. 


_  (1) 


We  will  now  examine  the  case  of  testing  two  observed  sets  of  estimates, 

_  (2) 

and  x  ,  for  equality  of  mean  vectors,  i.e.,  HQ:  or. 


(1) 

P1 

(1) 

p2 

it 

(2) 

M1 

(2) 

u2 

.  * 

(In  the  self-correlation  algorithm,  superscript  (1)  quantities  may  be  thought 
of  as  the  New  Data  and  superscript  (2)  as  the  Rase  Data.l 
For  case  3,  to  test  Hq,  the  statistic 


T 


2 


¥2 

VN2 


_0)  _(2)  ,  _(1)  _(2) 

(x  -x  )'  T.  (x  -X  ) 


is  used,  where  1  =  ( o . ^ ) ;  =  ct|^  =  0^  by  assumption. 


This  statistic  follows  a  noncentral  distribution  with  p  degrees  of 
freedom  and  non-centrality  parameter 

M,  p  P 


6  = 


?5f-  ?  !  .M"  -  “!2X”  -j2,>- 

N1  N2  i=l  j=l  1J  1  1  J  J 
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For  case  4 


2  .0)  _(2)  _-l_(l)_(2) 

T  =  (x  -x  )'  I  (x  -x  )  is  used. 


where 


z  =  (o..)  and  o..  *  i  a..  ^  ^  o..  This  statistic  follows 

1j  1J  N.  N,  1J 


a  non-central  x  distribution  with  p  degrees  of  freedom,  and 
non-centrality  parameter 

P  P  (D  (2)  (1)  (2) 

6  s  £  £  .  (is  -Ui  )(y,-  ) 

i  =  l  j  =  l  U  1  1  J  J 

The  preceding  results  (cases  3  and  4)  can  be  found  in  Johnson  and  Leone 
(1964, pp. 296-297). 

Case  5  is  analogous  to  case  3  except  that  I  must  be  replaced  with 


S  = 


Then  T  becomes 


N.N 


T2.  -ICiL 


N,  *  n2 


Appendix  IV  derives  the  distribution  of  this  statistic,  which  leads  to  the 
critical  region 

(N1  +  N2-2)p 


t2> 


Nj  +  N2  -p-1  P,N1  +  N2  'P_1 


(.0 


where  a  is  the  level  of  significance.  It  is  only  in  case  6,  unequal  and  unknown 
covariance  matrices,  that  there  is  not  a  precise  statistical  procedure 
yet  developed  to  test  the  equality  of  mean  vectors.  In  this  case,  other  methods 
such  as  computer-intensive  data  handling  techniques  must  be  used. 
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The  preceding  tests  Mill  yield  precise  results  only  when  all  of  the 
assumptions  concerning  the  data  are  met.  However,  in  some  applications  not 
all  of  the  data  will  be  normally  distributed  observations.  The  remainder  of 
this  report  discusses  techniques  used  to  analyze  the  validity  of  the  tests  when 
assumptions  are  not  valid. 


Chapter  III  -  Multivariate  Distributions  with  Component  Dependence 

This  chapter  describes  the  mathematical  methods  used  to  generate  data 
to  test  the  robustness  of  the  statistical  tests  of  Chapter  II.  The  robustness 
results  are  contained  in  Chapter  VII  while  the  computer  techniques  used  are 
in  Chapter  VI. 

To  simulate  relevant  data,  it  is  important  to  understand  the  characteristics 
of  the  observations  made  by  a  sensor.  Examples  of  these  observations  are  the 
signal  parameter  and  location  vectors.  A  signal  parameter  vector  may  have 
four  characteristics  or  components:  pulse  repetition  interval,  pulse  width, 
scan  rate,  and  frequency.  The  location  vector,  on  the  other  hand  ,  may  have 
longitude,  latitude,  and  altitude  as  its  components.  Each  data  vector 
is  independent  of  the  others  since  the  sensor  takes  readings  from  different 
locations  at  different  times;  however  due  to  the  nature  of  the  type  of  data 
collected,  the  components  of  each  data  vector  may  not  be  independent  of  each 
other.  A  good  example  of  such  dependence  in  the  signal  parameter  vector  is  the 
relationship  between  the  pulse  width  and  the  pulse  repetition  interval.  These 
two  components  are  related  through  the  peak-to-average  power  ratio. 

When  simulating  data,  we  assume  p=4.  That  is,  we  assume  the  mean 
vectors  have  4  components  and  the  covariance  matrices  have  dimension 
4x4.  (JPL  has  indicated  this  would  be  sufficient  for  the  desired  appl ications. ) 
The  Clinic  team  generated  data  from  the  multivariate  normal,  gamma,  and  lognormal 
distributions.  The  choice  of  these  particular  distributions  was  based  on  two 
criteria.  First,  these  distributions  have  reasonable  statistical  properties 
and  hence  the  mathematical  development  was  possible.  Secondly  they  have  the 
.desired  properties  for  testing  robustness.  The  multivariate  gamma  and 
lognormal  provide  examples  of  non-normal  or  skewed  distributions  while  the 
multivariate  normal  is  used  as  a  "control"  distribution. 


The  general  approach  in  generating  the  desired  multivariate  random 
samples  is  as  follows.  We  first  specify  the  desired  component  distributions 
and  component  dependency  (i.e.,  the  covariance  matrix,  denoted  by  >:  =  (o^)). 
To  obtain  the  specified  dependency  the  components  of  the  vectors  in  the 
random  sample  are  assumed  to  be  expressible  as  a  combination  of  appropriately 
chosen  univariate  distributions,  we  used  the  notation 


x  = 

i 

I 

and  the  y^  are 
In  matrix  form. 


x>\ 

X  = 

*2 

*  A  • 

y2 

X- 

J 

y3 

X4/  *4  J 

Thus  the  problem  becomes  one  of  finding  a  coefficient  matrix  A  such  that  the 
multivariate  distribution  x  will  have  the  specified  distribution  and 
covariance  matrix.  The  Clinic  team  named  the  process  of  finding  the  coefficent 
matrix  A  "backsolving"  since  the  algorithm  starts  with  the  desired  result 
(the  covariance  matrix  ,  z  )  and  ends  with  the  means  (the  coefficient  matrix, 

A)  by  which  the  problem  is  solved.  Once  the  coefficient  matrix  A  has  been 
found,  the  univariate  random  variables 


XA 

x2  1 


where  x.  =  z  a. .y.  ,  i=  1 ,. . . 
1  j=l  3 


j 


independent  univariate  random  variables. 


y.j  are  generated  using  standard 


techniques  which,  together  with  A,  produce  the  desired  multivariate 
random  samples. 


Before  looking  at  the  specific  distributions,  several  conments  should 
be  made.  First  since  the  covariance  matrix  is  symmetric,  we  need  only  work 
with  ten  of  its  entries;  therefore  we  assume  that  the  coefficient  matrix  is 
lower  triangular.  Secondly  when  testing  for  robustness  with  respect  to  the 
normality  assumption  in  case  5,  we  need  to  have  =  Eg  .  Thus  we  need 
to  generate  multivariate  random  samples  with  equal  covariances  matrices.  To 
do  this  we  first  show  that  we  can  backsolve  for  both  the  normal  and  gamma 
cases  provided  certain  relationships  among  the  entries  of  e  are  satisfied 
(see  equations  3.9  ,  3.10  and  3.11  ).  Since  backsolving  for  the  lognormal 

case  appears  to  be  very  difficult  mathematically,  we  avoid  it  by  choosing  our 
coefficient  matrix  A  so  that  the  resulting  covariance  matrix  can  be  backsolved 

m 

for  the  normal  and  gamma  cases.  Finally  we  note  that  although  the  Clinic  team 
approached  backsolving  using  basic  algebraic  manipulations,  Cholesky  decomposition 
techniques  could  have  been  used.  We  now  develop  the  backsolving  methods  used 
for  generating  multivariate  normal,  gamma  and  lognormal  data. 


NORMAL 

Constructing  specified  dependence  between  normal  random  variables  was 
relatively  simple  due  to  the  fact  that  any  linear  combination  of  normally 
distributed  random  variables  is  also  normally  distributed.  By  adding  four 
univariate  random  variables,  y.,  1*1,..., 4,  which  are  N(0,1),  new  random 
variables  were  created  as  follows: 


xi =  anyi  *  V2  *  Vs  *  V4  *  v  1  * ' . 4 


where  =  0  for  i  <  j 
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4 


I 


Each  is  distributed  normally  with  mean 
E  Dc^  =  V 

i  2 

Var  [x<]  =  r  a<t,  ,  and 
1  k*l  u 


Cov  ^  aikajk 

In  matrix  form,  the  new  equations  look  like: 
x  ■  A  y  +  u 

or  equivalently. 


(3.1) 


(3.2) 


(*? 

dll 

0 

0 

0 

fy\\ 

'1 

x2 

a21 

a22 

0 

0 

y2 

U2 

*  ' ' ' '  \ 

X  = 

X3| 

S 

a31 

3  32 

a33 

0 

• 

y3 

+ 

;x4j 

d41 

a42 

a43 

a44 

• 

...  -  - 

Given  a  certain  covariance  matrix,  r  =  (o-jj',  the  objective  is  t:o 
compute  the  coefficient  matrix,  or  the  A  matrix. so  that  the  cova  'iance 
matrix  of  x  is  equal  to  I  . 

Using  equations  3.1  and  3.2  ,  the  a^j  terms  can  be  computed.  For 
example,  by  equation  3.2  , 


.  i 


°  41  =  a41  an  +  a42  ai2  +  a43  d13  +  3 44  3  14 


s  a  a 
41  11 
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since  a]2  *  a]3  =  a]4  =  0.  And  by  equation  3.1  , 


2  2 

0|  =  1  »  or 


an  =  °i 


Thus, 


*41 


_  a41 


The  rest  of  the  a.,  terms  are  computed  in  a  similar  way  involving  numerous 

•  J 

calculations,  algebraic  manipulations,  and  substitutions.  They  are 
all  =  °1 


.  °12 


‘21  a 


1 


.  /detT 


22 


1 3 


31  o 


32 


a1  a23  -  q12°13 
°1  /det  B_ 


■/ 


det  C 


33  '  J  det  B 


>14 


41 


42 


CT1  °24  ~  0 1 2° 1 4 
°1  v/diO 
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The  gamma  distribution  was  chosen  to  represent  a  skewed  distribution. 
Through  its  parameters,  n  and  X,  the  gamma  has  a  great  deal  of  versatility. 
The  density  function  of  a  gamma  distribution  is 

V*>  ■  ^rTHT'e  'Xy 

Its  mean  and  variance  are 
E  [yj  =  j  and 

Var  [y]  =  A  respectively. 

To  simulate  multivariate  garrana  random  samples  we  assume,  as  in  the  normal 
case,  that  the  coefficient  matrix  is  lower  triangular. 

Let  x]  =  a]]y1 

x2  =  a21yl  *  aZ2y2 

x3  =  a31yl  +  a32y2  +  a33y3 

X4  =  a41yl  +  a42‘V2  +  a43y3  +  d44y4 

where  the  y.  are  independent  gamma  random  variables  with  parameters 
(ni »xi )  »  1  =  T*2*3*4-  i-e-,yi  ~  G(ni,xi). 


In  order  for  the  x^  to  be  distributed  as  gamma,  certain  restrictions 
have  to  be  imposed  on  the  n^  and  x^.  These  restrictions  become  clear  when 
we  look  at  the  moment  generating  function.  The  moment  generating  function 
(MGF)  of  y.  is 


Similarly, 


v-('-¥r 


the  other  hand,  the  MGF  of  Xg  is 


Mx  (t)  =  E  j^et^a21yl  +  a22y2^| 


r  r  a21 ty1i  .  a22ty2 
=  E  e  •  E  e 


a21 


A  'n9 


(- 1) 


(by  independence) 


A , 


1  2 

The  MGF  of  x~  indicates  that  x0  will  be  distributed  as  a  G(n,+n0,—  +— 
2  2  1  ^  a21  a2i 

only  if 


a21  _  a22 


(3.3) 


By  a  similar  analysis,  two  more  restrictions  are  imposed  for  x  ^  and  x^ 


to  be  gamma.  They  are 


a31  _  a32  _  a33 
A1  A2  X3 


(3.4) 


a41  !i2  =  ^43  a44 

A1  A2  A3  A4 


(3.5) 


In  sunmary,  if  restrictions  3.3  ,  3.4  ,  and  3.5  hold,  x^  will  be 


distributed  as 


r;  r  t 


The  mean  of  is 


4  nk 

E  [xj  =  z  —  a.k 
1  k=1  Ak  1K 


and  the  covariance  matrix  is 


a..  =  cov[x .  ,x.] 


1J 


=  E 


=  E 


i  J" 


/ 4  4  nk  \  /  4  4nk 

rZ  a..y.  -  z  T-a.A  z  a..y  -  z  —a, 

Li-1  lk  k  k=l  Xk  7\k  =  l  Jk  k  k  =  l  k 


4 


(“ik(yk  -  kkJ{k^,aj 


jk  lyk  “  Ak 


n 


i^l^k^k  ,2 


since  E 


IV 


y.  -  Hi' 

xj 


k  ^  j 
k  =  j 


Accordingly,  the  a.,  terms  can  be  expressed  in  terms  of  o--, 

1  J  1  J 


ij 

and  ni .  After  tedious  algebraic  manipulations 


aij 

=  0  ,  for 

a, , 

°11 

11 

/  o]1n1 

a21 

.  °21X1 

/  onn1 
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.Xmt  m  I  rn  r. ■  fi  m,A  mJrn  k  •  hm,»  ml  m  k  m  V  fc.A  m 


■-.V  .  ^  V  ... 


Similarly,  from  equation 


LOGNORMAL 


The  lognormal  distribution  is  a  skewed  distribution  which  is  related  to  the 
normal  distribution.  Given  that  y  is  normally  distributed  with  mean  «. 

p 

and  variance  a  then  x  =  exp[y]  is  defined  as  a  lognormal  distribution. 

The  expected  value  of  x  is 

p 

E[xJ  =  exp[w]  •  exp[o  J 

while  the  variance  is 

Var[x]  =  exp[2y]*exp[o2]‘(exp[a2-l]). 

Since  x  *  exp[yj,  the  condition  0<x<»  must  hold. 


Let 


Thus 


=  exp  La1]y1  +  b1  J 

x2  =  exp  La11y]  +  a22y2  +  b2J  • 

x3  ‘  exp  tallyl  *  a22y2  *  W3  *  b3] 

x4  *  exp  [allyl  +  a22y2  +  a33y3  *  a44y4  *  b4]  . 

where  the  a^  and  b^  are  scalar  constants,  and  the  y  are  Independent  standard.  ; 
nomal  distributions,  i.e.  y,  ~  H(0.1).  Each  x,  is  iognomally  distributed  ' 

i  • 

S,"Ce  Cjf,  ajjyj  *  biJ  ,S  distr,bljted  "Penally  with  mean  b,  and  variance  ' 

Using  the  model  shown  above,  the  expected  value  of  x.  is 

ELx.J  =  exp[b .]  .  expL  l  a  *  /  2] 

k=l  kk 
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The  variances  of  the  x..  ,  which  are  the  main  diagonal  elements  of  the 

covariance  matrix,  are 

Vor[x.j]  =  exp[2bi  +  (  akk)]  |exp(^akk)-l] 


Thus 


°i j  =  exp 


b,  +  b.  +  1/21 
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\k*l 


kk 
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l-  z 
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kk 
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2\  ' 

exp^ 

1  dkk) 
kk=l  kk' 

Please  see  appendix  V  for  the  actual  equations  for  each  a... 

^  *3 

The  relationships  needed  between  the  entries  of  the  covariance  matrix  z  =  (- 
in  order  to  backsolve  for  the  normal  and  gamma  cases  are 


c12  _  °T3  .  °13 

°22  °23  °24 

^23  _  °24 

°33  °34 

°33°44  , 

— r >  1 

°43 


Appendix  VI  shows  these  relationships  will  always  hold  for  a  multivariate 
lognormal  random  variable  defined  as  above. 
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Chapter  IV  -  Skewness  for  Multivariate  Distributions 


IS  BLANK 


The  previously  discussed  distributions  (Chapter  III)  were  chosen  in 

2 

order  to  study  the  robustness  of  the  Hotelling's  T  -statistics.  The  Clinic 

team  decided  to  analyze  the  effect  of  relaxing  the  assumptions  of  normality 

and  equal  covariance  matrices.  The  relaxing  of  equal  covariance  matrices 

is  discussed  in  Chapter  V.  To  analyze  robustness  with  respect  to  normality, 

the  normal  distribution  is  used  as  a  control  since  it  is  the  distribution 
2 

on  which  the  T  -statistics  are  based.  Two  descriptive  statistics  commonly 

used  in  the  literature  to  measure  the  non-normality  of  multivariate 

distributions  are  skewness  and  kurtosis.  Skewness  is  a  measure  of  how 

symmetric  a  distribution  is.  A  symmetric  distribution,  such  as  normal,  will 

have  a  skewness  of  zero.  Kurtosis,  or  excess,  is  a  measure  of  the  probability 

density  in  the  tails  of  a  distribution.  The  normal  distribution  is  sa  ^  to 

have  a  standard  level  of  kurtosis.  By  using  these  two  measures  in  statistical 

tests,  one  can  determine  whether  a  distribution  is  normal  or  not  [Kres  (1983, 

Table  26)]  These  tests  for  multivariate  normality,  developed  by  K.  V.  Mardia, 

were  not  studied  by  the  Clinic  team;  however,  Mardia's  theories  of  multivariate 

skewness  and  kurtosis  were  used  extensively.  In  a  Monte-Carlo  study  of  the 
2 

T  -statistic  in  a  univariate  case,  Mardia  found  that  the  level  of  skewness 
had  a  much  more  significant  effect  on  the  test  statistic  than  the  level 
of  kurtosis  [Mardia  ( 1 975 ,p. 1 67 ) J .  Similar  results  are  indicated  throughout 
the  literature. 

The  following  example  will  Illustrate  the  effects  of  skewness.  Figure  8 
shows  a  non-normal  skewed  distribution  centered  at  the  correct  meah  and  the 
assumed  normal  distribution.  This  case  can  lead  to  two  types  of  errors. 

If  a  candidate  entity  falls  outside  the  assumed  acceptance  region,  it  is 
considered  a  new  entity  and  is  added  to  the  Base  Data.  Thus  when  the  candidate 


Actual  Acceptance  Region 


Assumed  Acceptance  Region 


Skewed  Distributions 
Figure  8 


entity  falls  in  Region  I,  an  error  will  be  made  since  the  assumed  acceptance 
region  is  too  small.  This  error  will  create  a  phantom  entity  in  the  Data  Base 

The  opposite  problem  accurs  in  Region  II.  Since  the  assumed  acceptance 
region  is  too  large  compared  to  the  actual  acceptance  region,  the  tests  will 
incorrectly  associate  a  candidate  entity  with  a  known  entity  in  the  Data  Base. 
Hence  a  new  entity  will  not  be  detected  when  it  enters  the  battlefield. 

To  analyze  the  effect  of  skewness,  the  gamma  and  lognormal  distributions 
were  examined  along  with  the  normal.  The  garrma  and  lognormal  are  both  skewed 
distributions  common  to  statistical  analysis.  This  chapter  will  develop 
the  skewness  model  in  general  terms  and  then  present  the  skewness  formulas 
for  the  gamma  and  lognormal  distributions  used  in  the  Clinic  teams's 
analysis. 

Most  distributions  can  be  characterized  by  their  moments  about  the 
mean.  Specifically,  the  third  moment  is  the  measure  of  skewness: 
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^[(x-±l)'E  1  (y-u ) 


(4.1) 


Here,  x  and  y  are  independently  and  identically  distributed  random 
variables,  and  p  is  the  number  of  characteristics  in  each  vector.  For  p=2, 
the  following  theorem  gives  an  alternate  expression  for  6 
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Let 


a  =  ((xl~yl^°  +  (x2~ )(yi~vii )  and 


b  =((xrVJl)a12  +  (x2-ji2)a22)(y2-W2) 
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The  elements  of  E-1  can  be  expressed  in  terms  of  p  .  They  are 
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And  by  equation  4.3, 
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By  equations  4.3,  4.5  -  4.7, 
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The  same  procedures  also  apply  in  calculatinj  E[a  b].  First 
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Its  expected  value  is 
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It  can  be  expressed  in  terms  of  p  such  that 
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The  final  substitution  yields 
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Similarly, 
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and  therefore 


1  “  2  2 
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Finally,  results  from  4.8  ,  4.S  ,  4.10  ,  and  4.11  are  substituted  hack  into 

4.4  to  give 
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Regrouping,  we  obtain  the  desired  result: 
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This  result  corresponds  to  the  work  of  Mardia  in  the  two  dimensional 
case.  [Mardia  (1970,  p.  523. X\  Mardia  gives  a  general  formula  for  the  computing 
of  theoretical  skewness  of  a  multivariate  random  variable.  To  check  Mardia's 
results,  the  Clinic  team  developed  a  formula  for  the  case  P  =  4,  This  derivation 
is  presented  in  Appendix  VII.  The  team  then  implemented  both  Mardia's  formula 
and  the  team's  in  a  computer  program  to  compare  results  for  p  =  4.  It  was 
found  that  Mardia's  general  skewness  formula  is  correct.  The  formula  is: 
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rst  Mr's,t' 


Thus,  for  a  random  sample  from  a  multivariate  population,  S  or  1  is  computed, 
and  the  theoretical  moments  inherent  in  each  specific  distribution  are  used 
to  aet  a  measure  of  theoretical  skewness.  It  is  clear  that  if  the  distribution 
type  changes,  so  will  the  moments,  and  hence  the  value  of  skewness.  The 
theoretical  measures  of  skewness,  including  the  derivation  of  moments,  for 

the  multivariate  normal,  gamma,  and  lognormal  distributions  constructed  in 
Chapter  III  are  presented  next. 


NORMAL 


The  normal  distribution  is  symmetric,  and  has  a  third  moment  equal  to  'a 

zero;  thus,  its  measure  of  skewness  should  also  equal  zero. 

Given: 

X1  =  allyl  +  »*1 

x2  =  a21yl  +  a22y2  +  w2 

x3  =  a31yl  +  a32y2  +  a33y  3  +  M3  « 

x4  =  a41yl  +  a42y2  +  a43y3  +  a44y4  +  y4 


It  can  easily  be  shown  that  the  following  elements  of  Mardia's  skewness  equation  1 
which  apply  to  the  multivariate  normal  retribution  are  all  equal  to  zero. 
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The  first  and  third  terms  of  the  equation  on  the  right  are  equal  to 
zero  since  the  mean  and  third  moment  of  y..  are  equal  to  zero.  Due  to 
independence  among  the  y^  ,  and  that  E  [y^]  =  0  ,  the  second  term  will 
also  be  equal  to  zero  . 


In  equation  4.13, 
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Once  again,  every  term  on  the  right  side  of  the  equality  sign  is  equal  to 
zero  since  the  mean  and  the  third  moment  of  each  yi  are  equal  to  zero, 
and  all  the  y's  are  independent.  Similarly, 


E  [(x.  -  Vj)(xk  -  uk)j  =  0  . 

Substituting  equations  4.13,  4.14,  and  4.15  into  the  formula  derived  for 
01,4  in  APPendix  VII,  it  can  easily  be  seen  that  e,  =  0  for  p  =  4. 
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The  third  moment  of  x.  is 
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Since  the  are  independent,  the  interaction  terms  are  equal  to  zero 


(E  [w.]  =  0).  Thus, 
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Similarly, 
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and  for  i  <  j  <  k. 
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LOGNORMAL 


Given  that  y  is  n(0,<j  )  ,  then  x  =  exp[y]  is  a  lognormal  distribution 
and  the  coefficient  of  skewness  is 


E  L(x-T.E(x)r] 

(var[xj3^) 


For  computational  purposes,  the  important  elements  of  Mardia's  skewness  equation, 
e-j  are  based  on  the  following  dependencies  among  the  components  of  the  * 

x  vector: 

Xj  =  exp[a]1y1  +  b^  ‘  « 

x 2  =  exp[any1  ♦  a2^y2  +  b2] 

x3  =  exp^allyl  +  a22*2  +  a33y3  +  b2J 

X4  -  *  a22y2  *  a33y3  *  a^  .  b,] 

rfhere  the  y..  are  independently  identically  distributed  N(0,1).  -  — 

The  elements  of  Mardia's  skewness  equation 
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These  measures  of  skewness  are  specific  to  each  distribution  type 

since  each  distribution  has  different  values  for  the  moments  needed  in  the 

skewness  equation.  These  results,  then,  are  theoretical.  It  is  also 

possible  to  obtain  an  indication  of  sample  skewness,  b,  ,  which  is 

*  »P 

not  dependent  on  distribution  type.  This  measure  is:  • 

N  N 

b,  p  *  -4  E  £  L(xi  -  x)'S_1(x.  -  x)]3 
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.  • 

b1  can  be  tested  against  zero  to  see  if  the  sample  is  from  a  symmetric 

distribution.  This  test  is  presented  in  [Kres,  (1983, Table  26 ) J .  By  using 

both  the  sample  and  theoretical  measures  of  skewness,  the  Clinic  team  examined 

2 

the  robustness  of  Hotelling's  T  statistic  with  respect  to  distribution 
type. 
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Chapter  V  -  Unequal  Covariance  Matrices 

Another  assumption  that  the  team  relaxed  was  that  of  the  equality  of 

the  covariance  matrices.  It  is  important  to  analyze  the  robustness  of  the 
2 

T  -statistic  when  this  assumption  fails  because  there  is  no  statistic  derived 
for  the  case  of  unknown  and  unequal  covariance  matrices  (Case  6).  If  the  test 
statistic  for  case  5  is  robust,  then  it  would  be  acceptable  for  use  in  case  6; 
if  the  statistic  is  not  robust,  however,  using  the  wrong  statistic  would 
lead  to  phantoms  or  other  types  of  misinformation  in  the  Base  Data.  To 
illustrate  this  situation  in  the  univariate  case.  Figure  9  shows  two  symmetric 
distributions  with  the  same  means,  but  the  actual  distribution  of  the  data  has  a 
much  greater  variance  than  is  assumed,  i.e.  the  actual  distribution  is  "fat¬ 
tailed."  The  effect  of  this  incorrect  assumption  is  an  acceptance  region  which 
is  too  small  compared  to  that  of  the  actual  distribution.  As  in  the  case  of 
Region  I,  in  Figure  3,  phantoms  may  occur  because  an  estimate  that  falls  in  the 
shaded  region  is  from  a  known  entity,  but  the  test  will  indicate  that  any 
estimate  that  falls  outside  the  assigned  acceptance  region  is  a  new  entity. 

The  Clinic  team  altered  the  equality  of  covariance  matrices  by  allowing 
the  input  of  an  additional  positive  definite  matrix  to  be  added  to  the  covariance 
matrix  of  the  test  distribution.  In  this  way  it  is  possible  to  make  either  a 
larger  or  a  smaller  covariance  matrix  for  the  test. 

Chapter  VII  contains  the  results  obtained  by  relaxing  the  assumption 
of  equal  covariance  matrices. 


Chapter  VI  -  Computer  Simulation 


The  programs  written  for  the  Claremont  Graduate  School  Mathematics  Clinic 

under  the  auspice^ of  the  Jet  Propulsion  Laboratory  are  designed  to  simulate 

data  from  various  probability  distributions  and  perform  the  Hotelling's  T  -tests 

to  determine  whether  or  not  two  sets  of  data  are  from  the  same  distribution. 

2 

As  described  in  the  report,  the  T  -test  assumes  that  the  data  from  the  two 
populations  are  normally  distributed.  This  program  is  used  to  study  the 
robustness  of  the  test  when  the  data  is  not  properly  distributed;  in  particular, 
when  the  data  follows  a  skewed  distribution,  i.e.  a  gamma  or  lognormal.  The 
mathematics  of  these  tests  have  already  been  presented.  This  chapter  constitutes 
a  "user's  manual"  for  the  programs  written  by  the  Clinic  team. 

Since  simulation  requires  a  great  de*l  of  computer  time,  the  program 

N 

INPUT  was  written  to  create  input  data  for  the  simulation  program,  SIMULATE. 

With  data  generated  from  INPUT,  SIMULATE  is  not  interactive  and  requires  no 
user  attention. 

Specifically,  INPUT  sets  all  of  the  necessary  parameters  for  SIMULATE  to 
run.  These  parameters  are  as  follows: 

1.  The  output  file  for  the  information  from  the  simulations. 

2.  The  distribution  types  for  the  test  distribution.  Here  we 

have  l=Uniform;  *2=Normal;  3=Exponential ;  *4=Gamma;  *5=Lognormal ; 
6=Weibull;  *7=Cauchy;  *8=Gamma  and  Lognormal  with  the  same 
covariance  matrix.  The  starred  numbers  are  implemented  completely 
in  SIMULATE.  The  other  distributions  are  not  complete  and  are 
left  for  further  development.  The  default  is  a  Normal  distribution, 
i.e.  if  no  distribution  is  specified  a  normal  distribution  will 
be  used. 
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3.  The  test  number  depending  on  the  assumptions  made  about  the 
covariance  matrices.  The  test  number  is  1  if  the  covariance 
matrices  are  known  and  equal,  2  if  the  covariance  matrix  is 
estimated  from  the  test  data,  and  5  if  the  covariance  matrices 
are  unknown  and  equal.  Tests  3,  4,  and  6  require  further 
work  and  are  not  implemented. 

4.  The  parameters  necessary  to  generate  random  observations  from 
the  specific  distribution.  The  specific  method  for  obtaining 
observations  will  be  discussed  later. 

5.  The  amount  by  which  the  theoretical  means  of  the  Normal  base 
distribution  should  be  varied  from  those  of  the  test  distribution. 

6.  The  positive  definite  matrix  by  which  the  covariance  matrix  of 
the  base  distribution  should  be  varied  from  the  test  distribution. 

7.  The  theoretical  mean  of  the  base  distribution. 

8.  The  theoretical  covariance  matrix  of  the  base  distribution. 

9.  The  confidence  level  (1-significance  level)  for  the  test  where 
1 =99%  and  2=95%. 

10.  The  sample  size  of  the  generated  distributions.  A  large  sample  of 
100  or  a  small  sample  of  20  are  the  sizes  used.  If  other  sample 
sizes  are  desired,  then  it  is  necessary  to  alter  the  critical  test 
values  in  the  code  to  reflect  the  sample  size. 

11.  The  number  of  simulations  using  the  parameters  as  described  in 
2  through  10  above. 

INPUT  allows  the  user  to  input  data  for  as  many  runs  as  desired  so  SIMULATE  may 
generate  many  sets  of  results.  Given  the  above  parameters,  SIMULATE  may  generate 
the  prescribed  data  and  perform  the  specified  tests. 

The  flow  of  control  of  SIMULATE  is  straightforward.  The  parameters  are 


read  in,  the  observations  from  the  specified  distributions  are  generated,  and 
the  particular  test  is  performed.  Clearly,  there  are  many  details  to  be 
filled  in,  but  rather  than  describing  every  line  of  code,  a  description 
of  each  procedure  is  given  here. 

PROCEDURE  RANDOMIZE 

This  procedure  generates  a  seed  integer  to  be  used  subsequently  in 
the  function,  RND,  described  below.  This  seed  is  generated  by  calling  the 
external  system  Math  Library  function  FORSSECONDS.  The  output  of  FORSSECONDS 
is  a  real  value  corresponding  to  the  current  time  given  by  the  computer's 
internal  clock.  Since  this  value  is  given  in  milliseconds,  it  is  improbable 
that  the  same  seed  will  be  returned  twice.  This  greatly  improves  the  "randomness" 
of  this  random  number  generator. 

FUNCTION  RND 

This  function  is  called  throughout  the  program  to  return  a  random  number 
in  the  range  (0,1).  After  every  500  calls,  RND  will  generate  a  new  seed  to 
be  used  in  subsequent  calls.  RND  calls  the  external  system  Math  Library 
function  MTH$RANDOM  to  return  the  random  number.  See  Appendix  IX. 

PROCEDURE  INPUT  PARAMS 

INPUT_PARAMS  is  the  procedure  which  gets  all  of  the  necessary  information 
from  the  file  created  by  INPUT.  The  actual  parameters  read  in  differ  according 
to  the  distributions  to  be  generated. 

PROCEDURE  MAT_0UT_0NE 

This  procedure  outputs  an  NxN  matrix  to  the  output  file. 

PR0CEDUREMAT_0UT_TW0 

This  procedure  outputs  two  NxN  matrices  to  the  output  file,  side  by  side. 
PROCEDURE  INVERT 

This  is  an  external  procedure  written  in  BASIC  which  computes  the  inverse 
of  the  matrix  passed  to  it.  The  BASIC  matrix  handling  capabilities  are  much 


easier  to  use  than  trying  to  invert  a  matrix  using  brute  force.  At  this 
point,  no  investigation  has  been  made  concerning  the  accuracy  of  the  BASIC 
inverse  routine.  If  a  better  routine  is  found  (such  as  an  IMSL  routine) 
it  could  replace  the  current  inversion  procedure. 

FUNCTION  MAT  MULT  ROW  COLUfIN 

A  row  vector  and  a  column  vector  are  sent  to  this  function,  and  the 
scalar  result  is  returned. 

FUNCTION  MAT_MULT_ROW_MA.TR IX 

A  row  vector  and  a  matrix  are  sent  to  this  function.  The  result,  a 
row  vector,  is  returned. 

FUNCTION  TRANSPOSE 

This  function  returns  the  row  vector  which  is  the  transpose  of  the  column 
vector  passed. 

FUNCTION  DETERM 

This  function  computes  the  determinant  of  a  4x4  matrix  along  with  its 
sub-determinants.  The  particular  determinant  calculated  depends  on  whether  a 
2,  3,  or  4  is  passed  into  the  parameter  SIZE.  Each  calculation  is  done  by  follow 
the  general  formula  for  the  determinant  of  a  matrix.  No  recursion  is  used. 
PROCEDURE  UNISETS 

This  procedure  generates  observations  from  a  bivariate  Uniform 
distribution  with  parameters  (A1 ,  Bl)  and  A2,  B2).  This  procedure  was  inserted 
at  the  beginning  of  the  program  development  and  was  recently  updated  to  build 
4-characteristic  samples. 

PROCEDURE  NORMSETS 

This  procedure  generates  observations  from  a  Normal  distribution  with  tour 
characteristics.  This  procedure  requires  the  mean  about  which  the  observations 


are  centered  and  the  matrix  which  describes  the  dependencies  between  the 
characteristics. 

PROCEDURE  EXPOSETS 

Like  UNISETS,  this  procedure  is  incomplete.  It  generates  independent 
observations  from  an  exponential  distribution  with  parameter  x. 

PROCEDURE  GAMMASETS 

This  procedure  generates  4-variate  observations  from  a  Gamma  distribution 
using  a  sum  of  exponential  random  variables. 

PROCEDURE  LOGNORMAL 

This  procedure  generates  4-variate  observations  from  a  lognormal  distribution. 
The  methodology  used  to  build  dependence  between  the  parameters  described  in 
Chapter  III  is  used  here.  The  local  procedure  AMAKER  is  used  to  generate  an 
A  matrix  and  a  B  vector.  These  entities  are  generated  randomly  between 
specified  upper  and  lower  bounds.  These  are  used  to  generate  the  covariance 
matrix  which  must  be  positive  definite  in  order  to  solve  for  the  dependence 
matrix  required  to  generate  the  data  from  the  normal  distribution  used  as  the 
base.  A's  and  B's  are  generated  until  the  covariance  is  positive  definite. 
PROCEDURE  CAUCHYSETS 

This  procedure  generates  4-variate  observations  from  a  Cauchy  distribution. 

This  distribution  has  "fat-tails"  and  may  be  used  to  test  the  robustness  of 
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the  T  -squared  test  on  a  fat-tailed  distribution. 

PROCEDURE  WEIBULL 

This  procedure  may  be  used  to  generate  data  from  a  Weibull  distribution, 
but  it  is  not  completely  coded. 

PROCEDURE  STAT  TESTS 

This  is  the  main  procedure  of  the  program.  It  performs  the  statistical 
tests  described  earlier.  The  local  procedures  GAMSKEW  and  LOGSKEW  compute 


the  theoretical  skewness  using  the  moments  of  the  distributions  according 
to  Mardia.  Within  STAT_TESTS,  a  count  of  the  number  of  acceptances  and 
rejections  is  kept  to  summarize  the  results  of  the  simulations.  This  procedure 
also  computes  the  estimated  skewness  and  kurtosis  of  the  test  distribution. 
PROCEDURE  INFOJDUT 

This  is  the  output  procedure  of  the  program.  It  outputs  all  of  the  vital 
information  needed  to  analyze  the  test  statistics.  This  information 
includes  the  test  type,  the  test  distribution,  the  mean  vector,  and  the 
covariance  matrices  of  the  test  and  base  distribution,  the  parameters  used 
to  generate  the  data,  and  the  dependence  matrix.  After  the  summary  information 
is  printed  once,  the  user  can  suppress  repeating  the  summary;  in  this  case, 
only  information  that  changes  is  reported  after  the  summary. 

As  mentioned  above,  the  flow  of  SIMULATE  is  rather  straightforward.  The 
only  section  of  code  that  may  seem  unclear  is  case  8  in  the  FOR  Z:=l  TO  RUNS  loop. 
When  the  program  was  originally  developed,  it  was  designed  to  generate  data  for 
only  one  test  distribution.  When  the  decision  was  made  to  run  LOGNORMAL  and 
GAMMA  data  with  the  same  covariance  matrix,  it  was  easier  to  simply  generate 
lognormal  date,  perform  the  tests,  generate  the  gamma  data,  and  then  let  the 
program  continue  normally.  Thus,  the  necessary  code  was  exactly  duplicated  for 
running  the  lognormal  tests.  This  is  not  necessarily  efficient,  but  it  was  the 
best  solution  given  the  time  constraints. 

As  a  final  note,  the  reader  may  refer  to  Appendix  VIII  for  the 

simulation  techniques  for  generating  the  univariate  random  samples. 


Chapter  VII  -  Robustness  Results 

The  self-correlation  stage  uses  statistical  tests  to  check  if  one  set 
of  data  (New  Data)  is  from  the  same  population  as  another  set  of  data  (Base 
Data).  The  tests  are  based  on  assumptions  about  the  data,  such  as  distribution 
type  and  dependency  relationships.  The  study  of  robustness  concentrates  on 
a  statistical  test's  ability  to  perform  as  expected  even  though  the  data  does 
not  conform  to  the  assumptions  upon  which  the  test  is  based;  specifically,  if 
the  data  from  the  battlefield  are  not  normal,  or  the  assumptions  on  the 
covariance  matrices  are  incorrect,  will  the  tests  still  perform  accurately?  These 
assumptions  were  addressed  by  the  implementation  of  the  mathematical  techniques 
and  simulation  program  previously  described  (Chapters  III-VI). 

The  first  way  in  which  the  clinic  n  analyzed  the  robustness  of  Hotelling's 
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T  -statistics  was  through  the  use  of  non-normal  distributions.  As  mentioned  in 
Chapter  IV,  skewness  has  a  detrimental  effect  on  the  performance  of  these  tests 
in  the  univariate  case.  The  team  used  backsolving  techniques  (Chapter  III)  to 
generate  multivariate  skewed  distributions  to  be  used  to  test  the  T^-statistics 
for  cases  1,  2,  and  5.  (These  were  the  only  cases  simulated  since  case  3  and 
case  4  do  not  appear  to  have  relevent  applications  and  no  test  statistic  was 
derived  for  case  6.)  A  description  of  the  6  cases  is  contain  in  Figure  7, 

Chapter  II.  There  are  two  ways  that  the  tests  could  fail  to  perform  properly; 
the  test  could  reject  the  null  hypothesis  when  the  null  hypothesis  is  really 
true  (Type  I  error),  or  it  could  accept  the  null  hypothesis  when  it  is  really 
false  (Type  II  error). 

The  Clinic  team  tested  for  Type  I  error  by  simulating  Base  Data  and 
New  Data  with  equal  mean  vector  and  covariance  matrices.  At  a  5%  signigicance 
level  the  null  hypothesis  should  be  accepted  approximately  95%  of  the  time  when  it 
is  true.  Figure  10  shows  the  results  of  the  simulation.  Since  the  statistic. il  t.--. 
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Simulation  Results  1 
Figure  10 


are  based  on  the  assumption  that  the  data  follow  normal  distributions,  they  should 
be  accurate  when  normally  distributed  data  is  used.  That  is,  the  null  hypothesis 
should  be  accepted  approximately  95%  of  the  time  when  it  is  true.  The  first  row 
shows  that  in  each  of  the  three  cases  the  percentages  are  as  expected.  Since 
the  gamma  and  lognormal  distributions  are  skewed  distributions  (i.e.,  not  normal) 
we  would  expect  the  tests  to  fail.  The  second  and  third  rows  indicate  that  the 
tests  do  fail  for  skewed  distributions.  Except  for  the  lognormal  distribution 
in  case  5,  neither  the  gamma  nor  the  lognormal  distributions  lead  to  the 
acceptance  of  the  null  hypothesis  more  than  63%  of  the  time  when  it  is  true; 
in  fact  the  gamma  distribution  never  did  better  than  28%.  The  reasons  why  the 
lognormal  distribution  performed  well  for  case  5  are  unclear  at  this  point  in 
time.  It  is  suspected  that  the  use  of  pooled  estimates  for  covariance  matrices 
from  both  the  New  Data  and  the  Base  Data  could  have  compensated  for  the  effect 
of  a  slightly  skewed  distribution. 

The  Type  II  error,  accepting  the  null  hypothesis  when  it  is  not  true,  was 
tested  by  generating  data  with  equal  covariance  matrices  but  unequal  mean 
vectors.  Hence  we  would  expect  the  null  hypothesis  of  equality  of  mean 
vectors  to  be  rejected  a  large  percentage  of  the  time.  (That  is,  the  null 
hypothesis  would  be  accepted  a  small  percentage  of  the  time.)  As  shown  in 
Figure  11,  the  null  hypothesis  was  accepted  0%  of  the  time  for  all  three 
distribution  types  in  cases  1,  2,  and  5.  Thus,  for  Type  II  error  simulation, 
the  statistical  tests  performed  well  for  both  the  normal  and  skewed  distributions. 
In  other  words,  the  statistical  tests  appear  to  be  robust  with  respect  to 
Type  II  error  for  the  distributions  considered. 

The  final  simulation  performed  generated  data  from  the  three  distributions 
with  equal  mean  vectors  but  with  different  dependency  relationships;  thus 
pn  =  pb  but  t  ig.  This  simulation  was  used  to  check  the  robustness  of 
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equal  covariance  matrices  in  case  5.  Once  again,  with  equal  mean  vectors 

we  would  expect  the  null  hypothesis  to  be  accepted  approximately  95%  of  the 

time.  Figure  12  shows  the  results  of  this  simulation.  While  the  normally 

distributed  data  performed  well,  the  two  skewed  distributions  did  not  (98% 

vs.  23%  and  71%).  We  note  that  these  results  are  similar  to  those  shown 

in  Figure  10  for  case  5.  Thus,  these  preliminary  results  seem  to  indicate 

that  small  to  moderate  differences  in  the  two  covariance  matrices  have  little 

effect  on  the  performance  of  the  tests.  These  results  imply  that  it  is 

the  skewness  of  the  distributions  that  causes  the  non-robustness  of  the  tests. 

We  note  that  it  is  important  to  analyse  robustness  of  equal  covariance 

matrices  since  there  is  no  test  statistic  derived  for  the  case  of  unknown 

and  unequal  covariance  matrices  (case  6).  If  the  test  statistic  for  case 

5  is  robust,  then  it  would  be  acceptable  for  use  in  case  6. 

In  summary,  the  Clinic  team  used  a  simulation  approach  to  test  the 

2 

robustness  of  Hotelling's  T  -statistics  with  respect  to  normality  and  equal 

covariance  matrices.  Three  separate  simulations  were  performed  involving 

normal,  garrma,  and  lognormal  distributions.  The  results  of  the  simulations 
2 

show  that  the  T  -statistical  tests  are  not  robust  with  respect  to  the  normality 
(or  symmetry)  assumptions.  In  other  words,  if  the  data  being  analyzed  comes 

O 

from  a  skewed  distribution  and  a  variation  of  Hotellings  T  -statistic  is 

used,  the  decision  made  by  the  test  will  probably  be  incorrect.  More  studies 

7 

need  to  be  made  on  the  relationship  between  skewness  and  the  T  -statistics. 

2 

Using  T  -statistics  that  assume  equal  covariance  matrices  on  datasets  that 
actually  have  slightly  different  covariance  matrices,  however,  seems  to 
have  an  insignificant  effect  on  the  performance  of  the  tests. 
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Figure  12 


Chapter  VIII  -  Concluding  Remarks 


In  this  report  the  Clinic  team  first  developed  the  mathematical  techniques 
and  simulation  program  required  to  generate  multivariate  distributions  with 
component  dependency.  These  results  were  then  used  to  investigate  the  effects 
of  skewed  distributions  and  unequal  covariance  matrices  on  the  statistical  tests 
in  the  self-correlation  algorithm.  (The  results  of  these  investigations  are 
contained  in  Chapter  VII.)  The  Clinic  team  plans  to  continue  working  in  this 
area.  In  particular,  since  preliminary  results  indicate  that  skewed  distributions 
cause  inaccurate  conclusions  in  the  statistical  tests,  the  relationship  between 
the  robustness  of  these  tests  with  respect  to  multivariate  skewness  and  kurtosis 
will  be  studied  further.  This  will  be  the  main  objective  of  the  next  report. 
Applications  of  Correlation  Techniques  Battlefield  Identification  I_^. 

There  are  at  least  two  other  areas  that  clearly  should  be  studied  in 
future  reports.  One  topic  would  be  the  investigation  of  the  robustness  of  the 
statistical  tests  with  respect  to  other  assumptions.  These  assumptions  include 

1.  The  sensor  data  is  unbiased. 

2.  The  error  term  of  the  mean  vector  is  independent  of  time. 

3.  The  New  Data  describes  only  one  entity. 

4.  The  entities  are  stationary. 

The  Clinic  also  plans  to  develop  a  final  computer  package  to  implement  results. 
This  package  would  be  a  user-friendly  process  which,  among  other  things, 
performs  goodness-of-fit  tests  on  the  data  to  determine  whether  or  not  the 
proposed  statistical  tests  are  appropriate. 
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Finding  the  Maximum  Likelihood  Function  for  y/z 


The  likelihood  function  L(£';r  )  is  given  by 

i  *  n 
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L  = 


«P  t-%  r  U  -  £)-  i''(x  -  ,)]  (AI . I ) 
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To  find  the  maximum  likelihood  estimates  for  ^  and  E  in  Equation  AI.l 


we  begin  by  taking  the  natural  logarithm  of  Equation  AI.l.  Denoting  £*  and  t* 


-1 


to  be  the  maximum  likelihood  estimates  for  y  and  E  respectively,  we  get 
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In  L  =  -%pN  ln(2ir)+55N  ln|y>*|  E  (x  -  y*)ij/  *(x  -  u*)  (AI.2) 

a=l  a  1 


where  y*  and  maximize  In  L. 

The  following  lemma  will  be  useful  in  solving  the  above  equation  for 
and  i^*. 


Lemma  1 :  Let  x  1 . x^  be  N  (p-characteristic)  vectors.  Then  for  any 


vector  b, 
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The  second  and  third  terms  are  equal  to  zero  since  z(x  -x)  =  Ex  jN( — r2— )  =  0. 


Thus,  let  b=y*  and  apply  Lemma  1,  to  obtain  that 
N  N  _ 
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where  A  =  r(x  -  x)(x  -  x)'. 
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Recall  the  properties  of  the  trace  of  a  matrix: 
If  C  is  mxn  and  0  is  nxm  then 


m  n 

tr(CD)  =  l  l  C^d..  =  tr(DC) 
i=l  j=l  J1 

tr(C+0)  =  tr(C)  +  tr(D) 

Applying  Equation  AI.3  and  the  properties  of  the  trace  of  a  matrix  to  the 

third  term  of  the  right  side  of  Equation  AI.2,  we  have 

N  N 

Z  (x  -**)'** (x  tx*)  -  tr  I 


11 

*  tr  Z 

a=l 

N 

-  tr  \|r*  Z  (2ta“J±*)(*a1i*)' 
a*l 

=  tr{\|r*[A+N(x-]i*)  (x-]i*)']  } 

*  tr[^*A]+tr  [\|r*N(x-ji*)  (x-ja*>  '1 
■  tr[\|r*A]+N(x-^i*) '^*(x~ti*) 

With  the  result  of  Equation  AI.4,  Equation  AI.2  can  be  rewritten  as 


(AI.4) 


In  L  -  -ljpNln(2TT)+^  ln|f*|-Htr(t*A)  -  ^(x-H*)  't*(x-p*)  (AI.5) 


We  only  need  to  maximize  the  second  and  third  terms  of  the  right  side  of 
Equation  (AI.5)  because  the  first  term  is  a  constant,  and  the  fourth  term  is 
equal  to  zero  then  p_*=x. 

To  maximize  the  second  and  third  terms,  which  are  '-sN  ln|^*|  and  btrU*A) 
respectively,  we  must  use  the  following  lemma. 
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wlie  re 


Lexma  2:  Let  f(c)  -*>N  ln|c|  -  h  >■  i  d 

i.j-1  U  Jl 

c  *  (c-A  is  positive  semidefinite  and  0  =  (d..)  is  positive  definite.  Then 

1  J  •  J 

the  maximum  of  f(c)  is  taken  at  C=ND'^  and  the  maximum  is 
f ( ND- 1 )  =  *spN (l n  N)-4N  In  |D|  -^pN. 

Applying  Lemma  2  to  Equation  AI.5  by  letting  C=i^*  and  D=A,  we  obtain 
In  L  =  -4pN  In  (2n)  +  bpN  In  N  -  In  |A|  -  4pN.  (Al.b) 

or,  when  the  second  and  third  term  of  the  right  side  are  combined,  we  obtain 

In  L  =  -!jpN  In  (2n)  +  JjN[ln  -yj^-  ]  -  bpN  (AT .  7) 

After  taking  the  exponential  of  Equation  AI.7,  the  result  is 

where  \hhH  =  ~hH  =  N*4pN  lAI"^  * 


APPENDIX  II 


2 

The  Distribution  of  T  for  Case  2 
Unknown  Variance-Covariance  Matrix,  H^:  psUq- 

2  _  -1  - 

We  wish  to  find  the  distribution  fo  T  =N(N-1  )(x-^)'A  (x-pg),  where 

x  is  distributed  as  N(±j0,|).  This  will  be  done  by  several  transformations 

o 

which  associate  T  with  a  simple  and  known  statistical  distribution. 

Remember  that  the  matrix  A  is  defined  as 


N 

A  ■  1  (x  -*) (x  -x)  ' 

_  —a  —  —a  — 

a**l 

N  /x  ,-x 

■  Z  )  (xal"X;’*,,,xairXp)  J"1, 

o*l  l  x  -x  '  K 

V  GP  Py 


-al-Y 


Note  that 


I  is  not  independent  because  each  x .  is  dependent  on  the 

J 


VX  -X  . 
ap  p' 


corresponding  x  ..  This  dependency  can  be  eliminated  by  summing  up  to  N-l, 
aJ  N-l 

i.e.,A  is  distributed  independently  as  E  z  z  z  is  independent  and  is 

ICt  Ct  Cl 

distributed  as  N(0,z).  We  denote  the  sample  covariance  matrix  by  S. 


o  _  I 

Our  first  transformation  is  to  let  T  -jr'S  where  y  =  /nTx-pq)  and  y 
is  distributed  as  N(x,z).  The  objective  of  this  transformation  is  to  simplify 
the  mean  so  that  it  will  become  zero  under  the  null  hypothesis.  In  notation. 


1  -  E[y] 


EIt/n(x-Uo)1  -i/iTEKx-^)] 

-V^(Elx]-p0) 


Var  1^1  *  Var  Q>  1  *  N  VaT  t  <2*-P-0>  3 

-  N(Var(x)+Var  (^)-Cov(x»iio^ 

»  N*Var(20*N*  -  £ 

Our  second  transformation  is  to  let  D  be  a  nonsingular  matrix  such  tlvt 
Di:D '= I ,  and  define 

2*  -  Dy 
S*  -  DSD' 

Y*  -  DY 

2  -1 
T  is  now  equal  to  y*'S*  y*, 

because  T  *  z's~  y 

-  (y 'l)S-1(iy) 

-  X'[D* (D*) -1 ]  S_1(D"1D)^ 

-  (Dy)'(DSD')-1(Dy) 

*  y*'S*  ^y* 


y*  is  distributed  as  N(y*,I). 

The  third  transformation  is  to  let  the  first  row  of  a  pxp  orthogonal 
matrix,  Q,  be  defined  by 


qli 


i"l * •  •  •  >  P 


In  other  words. 


‘1  * 


V  * 

y2 


\  Y*'Z* 


anything  so  that,  u  is  orthogonal 


j 


The  rest  of  the  matrix  Q,  other  than  row  1,  can  be  anything  as  long  as 


orthogonality  is  maintained. 


p  *2 
* 


*■  y 


>'•  qu  “1  because  — *7^ 


a=l 


1-1 


11 


P 

>;  y  *2 

i-1  1 


Q,  as  defined,  is  a  random  matrix.  What  we  want  to  do  is  to  express  T 
a  scalar  form  rather  than  in  a  matrix  form.  This  can  be  done  by  letting 


U  -  Qy* 

B  -  Q(N-1)S*  O' 


With  the  above  definitions,  the  first  component  of  U  now  becomes 

“i  ■  j^xA* 


Vii(,u^) 

1=1 


Jl*'l*  1  (qli) 
i=l 


•  1  *  v  i*'z* 


The  remaining  components  of  U  become 

p 

U.  «  L  q .  .v,*  j  i  1 


•  0 


Therefore,  k  *[  0 


p 


Then, 


because  U'B"1^  -  (Q^*)' (Q(N-i)s*Q')_1(c&*) 

-  (Q')”1]  ~~  [o‘1Q]z* 


-1 


(s*) 

Z  jjrj —  1*  (recall  y*-Dy  and  S*»DSD  ; 


(Dy) 


,  (DSD') 
N-l 


-1 


(Dy) 


y^D<(D<)~1S~LD~:LDv 

N-l 


.  r±h 

N-l 

Re-expressing  Equation  AII.l 

2 


T 

N-l 


1  1 


i^2  b11 


In  Anderson's  notation, 
1 


,-l 


11  “  bll  "  b(l)B22  b'(1) 

•  b„  2  „  where  B  A1 

11  2”"’P  \  C l)b!2, 


is  a  partitioned  matrix. 
Hence, 


2 


T 

N-l 


11 *2 , • # , p 
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(AII.l) 


2 

The  denominator  is  conditionally  distributed  as  x  (chi-square)  with  N-p 

degrees  of  freedom;  it  is  the  sum  from  1  to  N-p  of  the  square  of  w  where  w 

is  distributed  as  N (0 , 1 ) .  The  numerator,  on  the  other  hand,  is  distributed  as 
2 

a  noncentral  x  with  p  degrees  of  freedom  and  noncentrality  parameter 

2 

y* '  y *  =  Thus,  ^  ~  P  *- -  )  is  distributed  as  a  noncentral  F  wi 

p  and  N  -  p  +  1  degrees  of  freedom.  The  noncentrality  parameter  is  £  =  »  " 


APPENDIX  III 


Using  Likelihood  Ratio  Criteria  to 
Test  the  Equality  of  Variance-Covariance  Matrices 


Let  x  ^  be  an  observation  from  the  population  where 

— n 

ftsl....»Ng  and  g=T,...,q.  (q  is  the  total  number  of  populations.) 

x(9)  is  N(£^,xg)  snd  is  a  column  vector  of  size  p.  Let 

p  =  number  of  characteristics  , 

N  =  number  of  observations  in  the  gtn  population  , 

9 

N  =  total  number  of  observations  , 

A  -  (x  («)-7(»;))(x  (K)-i(«>>. 

Q  — n  —  —/I  — 


M  •  I 


vr  i , .  • .  ,<i  ,  and 


A  =  Sum  of  Ag's. 


We  want  to  test  Hq.-  *  •••  =  Xg 


The  likelihood  function  is 


q 

n 


g-l  (ZnJ^pNgl^l^g 


exp 


i-»s  z8  (x  (g)-1j.(8))r1(x  (8)-ix(g))] 

-a  c  g  -a 


a-1 


Define  n  as  the  parameter  space  where  each  xg  is  positive  definite  and 

is  any  vector.  Define  u  as  the  parameter  space  where  x-j  =  =  •••  =  Xc 

and  is  any  vector.  Then  the  maximizing  values  are 


;(9)  =  *  (g),£  ,  ja 


g  N 


over  ft 


and 


~(g)  B  i  (g)  :  .  a 


.  xg  =  over  * 


77 


Therefore,  the  likelihood  ratio  criterion  to  test  is 


q  A  N  /2 

n  -9-  9 

11  ki 


=  1 


q  *iN 

H  lAgl  9 

-  _ 

1AI* 


(AIII.l) 


5  OpN9 


x  <  A  (a)  where  \(a)  is  defined  such  that  Equation  AIII.l  holds  with  probabi  lity 
a  when  Hq  is  true. 
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APPENDIX  IV 

r2 


The  Distribution  of  T  for  Case  5 
Unknown  but  Equal  Variance-Covariance  Matrices 


H  .  CD  (2) 
H0  “J± 


This  derivation  uses  the  results  of  Appendix  II  to  show  that 


N„N 


12  ,-(1)  — (2) .  ,  -1,-(1)  -  (2) 


Nl+N2 


(xv  -xv S  (xv  -x 


) 


follows  an  F  distribution  with  p  and  N1  +  N2~p-1  degrees  of  freedom.  Recall 
that 


Set 


Nx+N2-2 


Zk  (x  (k)-x°°>(x 


k-1  a-1 


N1N2 


-  v  nl+n2 


(7(1)-x(2)) 


Under  the  null  hypothesis  z  is  distributed  as  N (0,i).  Defining  z 
this  way  accomplishes  the  same  goal  as  the  first  transformation  in  Appendix  II 
Therefore,  the  remaining  transformations  of  Appendix  II,  when  applied  to  z 
will  yield  the  desired  result.  Specifically: 


2  ,  <Ni+N2~2)p  (a) 


n1+n2-p- 


where  a  is  the  level  of  significance. 


APPENDIX  V 

Elements  of  Covariance  Matrix  of  a  Lognormal 


PREVIOUS  PAGE 
IS  BLANK 


*1  =  exp[a11y1  + 


*2  =  exPt^i'i'Xi  +  a  22^2  +  b2^ 

x3  =  exP^allyl  +  a22y2  +  a33y3  +  b3^ 

x4  =  exp[a1]y1  ♦  +  *33*3  +  a44y4  +  b4J 


and  the  y.  are  i.i.d.  N(0,1), 


i  i  2 

and  (  i  a-.yj  +  b.  ~  N(b,,  I  a-H  )  for  1-1, ..,4. 

j  =  l  JJ  J  1  1  j  =  1  JJ 

The  covariance  matrix  elements  are: 

°n  3  Lexp(2bi  +  a,,2)]  •  Lexp(an2)  -1] 

o22  =  [exp(2b2  +  a„2  +  a222)]  .  [expfa^2  +  a^2)  -1] 

°33  =  [exp(2b3  +  an2  +  a222  +  a332)]  •  [exp(an2  +  a^2  +  a^2)  - 


44  =  [exp(2b4  +  an2  +  a222  +  a332  +  a442)]  •  Lexp(an2  +  a222  +  a332  +  . 

2  2 

2an  +  a??  ? 

°12  =  °21  =  t-exp(bi  +  b2  +  — ^ - —  )J  •  [exp(a11^)  -1] 


’23  =  °32  =  ^rxP(D2  +  b3  + 


°?4  3  °4?  =  lexP(b9  +  *>A  + 


2a))2  *  Za222  *  a332 


)  ]  •  Lexp(a,,2  +  a222)  - 


2aH2  *  1&Z2  +  a33?  +  a342 


-)  J  •  [exp(an2  + 


[exp(b3  +  b4  + 
[exp(an2  +  a22 

[exp(b1  +  b3  + 


[exp(b1  +  b4  + 


APPENDIX  VI 

Covariance  Matrix  Restrictions  for  Lognormal  Case 

Proof  that  the  restrictions  on  the  covariance  matrix  also  hold  true  in  the 
lognormal  case .  Let 


Let 

'•u 

°1 2 

°1 3 

°14 

°21 

°22 

“23 

°24 

°31 

°32 

°33 

°34 

°41 

°42 

°43 

°44 

The  following  relationships  must  be  satisfied  in  order  to  backsolve  for 
the  Gamma  or  Normal  distributions. 


a33°44 


(AVI.l ) 


(AVI. 2) 


(AVI. 3) 


Using  the  results  of  Appendix  V,  AVI.l  becomes 

a  2 

[exp(b1  +  b2  +  a^2  +  -y-)  ]  *  [exp(a1]2)  -1] 

[exp(2b2  +  an2  +  a222)]  *  texp^all2  +  *22^ 

2  2  2 

Lexp(b1  +  b3  +  — ^ ^ - —  >3  *  [exp(an2)  -lj 

=  -  -  2 

2a, ,  +  2a?„  +  a.,.,  o  2 

[exp(b2  +  b3  +  — U— g - — - — )]  .  [exp(an^  +  a22£)  -1] 


[exp(b1  +  b„  + 


)J  *  Lexp(a11  )  -1] 


[exp(b2  +  b4  +  — '-!■ - 

“44 

— )]  •  [exp(an2  +  a2 

After  simplifying  each  ratio, 

2 

a  2 

a22 

exp  [b1  -  b2  -  2  ]  =  exp[b1  -  t>2  " 

a22 

2 

d 

]  =  exp  [b]  -  b2  -  - 

Clearly,  relationship  AVI.l  holds  in  lognormal 

case . 

For  the  relationship 

°23  °24 

°33  °34 

[AVI . 

2  2  2 
2a,/  +  2a  ‘  +  a  ‘ 

[exp  (b,  +  b,  +  — □ - 5-^ - ~ 

-  )J  • 

[exp(a,,  +  a„„  )  -1J 

[exp(2b3  +  a^ 


)]  *  [exp(a^  +  a22  +  a. 


[exp(b0  +  b„  + 


2  2 
a33  +  a44 


-  ]  •  [exp(a 
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For  the  relationship 

°33  °44 

- 2  >1,  (AVI. 3) 

°43 


^[exp(2b3  +  an2  +  a222  +  a332)j  [exp(an2  +  a222  +  a332)  -1  ]  • 
[exp(2b4  +  di  ] 2  +  a222  +  a332  +  a442J]  [expla^2  +  a222  +  a^2  +  a42: 


^[exp(b3  +  b4  +  an2  +  a222  +  a332  +  ~~  )]  •  [exp(an2  +  a222  +  a332) 

exp(an  +  a22  *  a33  +  a442)  -1  ^  ^ 

exp(a112  +  a222  ♦  a^2)  -1 


Since  all  three  relationships  are  valid,  backsolving  from  a  lognormally  generated 
covariance  matrix  to  a  normal  or  gamma  distribution, is  possible. 


APPENDIX  VII 
Derivation  of  s 


1,4 


Let  x  and  2  be  independently  and  identically  distributed  random 
variables  Define: 


by  matrix  multiplication. 


Let 


a  ■  w,  E  y.a 
1  i=l  1 


il 


b  =  w„  E  y.a 
c  i=l  1 


i2 


C  =  E  y.a 
J  i  =  l  1 


13 


d  =  w.  E  y.a 
*  i=l  1 


i4 


(AVII.l ) 

(Mil. 2) 

(MU. 3) 

(AVI  1 .4 ) 


The  problem  restated  is 

B1>4  «  E[(a  +  b  +  c  +  d)3]  ,  (AVI 1 .5) 

It  can  be  shown  that 

(a  +  b  +  c  +  d)3  =  (a3  +  b3  +  c3  +  d3)  +  (3a2b  +  3ab2  +  3c2d  +  3cd2 

+  3a2c  +  3ad2  +  3bc2  +  3bd2  +  3ac2  +  3a2d 

+  3b2c  +  3b2d)  +  (6abc  +  6abd  +  6av.d  +  6bcd)  (AVII.6) 


Note  that  equation  AVII.6  consists  of  three  groups  of  like  terms. 

Substituting  equations  AVII.l  through  AVII.4  into  AVII.6  and 
expanding: 


a3  = 


3  ^  i  1  3  4  4  i  1  ?  ■;  i  4  4 

w-,  [  E  (y.-a11)3  +  31  E  (yia1,)2y.aJ  +  6  E 

'  i=l  1  i=l  j=i  1  3  j  =  l  i=l 

j*i  U3 


1 1  i 

y^o  J 
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The  remaining  terms  in  the  second  group  of  equation  AVII.7  are  calculated 
similarly,  and 


???????????? 
3|Vb  +  a2c  +  a2d  +  tTa  +  b^c  +  b^d  +  c2a  +  c2b  +  c2d  +  da  +  d2b  +  d2c] 


=  3[  i  z  wi2w.(  e  *  z  *  1  Vky/lokl°tj)] 

i=l  j*i  1  J  £=1  k=l  k  1  £=1  m=l  k=l  rk  l 


4  4  4 


mi  ki  ej, 


Final ly, 


A  ii  4  i?  ^  k3 
abc  =  w,w„w  (  i  y.o  )(  z  y,oJ  )(  z  y.a  ) 
12  J  i=l  1  j=l  J  k=l  * 


(AVI  1 .8) 


444  i 1  12  k3 

=  w,w.,w,  z  z  i  y4y.ybo  oJa 

1  2  3  i=l  j=l  k=l  1  J  k 


Therefore,  the  last  group  is 

444  i 1  12  k3 

6Labc  +  abd  +  acd  +  bed]  =  6[w,w,w,  ill  y.:y.j.yi°  c  o 

1  2  Ji=l  j=l  k=1  1  J  K 

4  4  4  il  j2  k4 

+  w,w?w.  i  z  z  y,y  .y  a  o 

1  2  4i=l  j=l  k=l  1  J  k 


4  4  4  i  1  j  3  k4 

+  www  z  z  z  y.y.y.o  a  o 

1  J  4i=l  j=l  k=l  1  J  k 


4  4  4 

+  w^w.w.  z  z  z  y  .y  .yLoi2o^nk4  ] 


y  J  i-7  L  1 

4  4  4i=1  J=l  k=1  l  J  k 


(AVI  I . 9) 
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Substituting  equations  AVII.7,  AVII.8,  and  AVII.9  into  AVI  1.5,  and  using  the 
properties  of  the  expected  value  operator,  then 


4  3  f  4 

~  rf.  .  j  i  I 


E  E[w  3J  E  E[y.3](olk)3  +  3  E  E  E[y.Vj(o  K)  oJ 
=1  k  i  i=l  1  i=l  j=l  1  J 

J7i 


4  4 


2., 


Ik  2k  3k  ,  rr .  t  lk_2k  4k  A  cr .  i  Ik  3k  4k 


-  6  (ECy^lo' VVK  +  Etyy^o' +  Ety^V 


GO  + 


rr  ,  t  2k  3k  4k  v /i 
+  E[y?y3y4Jo  o  o  )}] 


,44  P  44  P  ki?H 
+  3  E  i  E[w.  w. J(  e  e  E[y,  y  ](o  )  oU 

li=l  j=l  1  J  2=1  k=l  K  1 
j^i 

4  4  4  _  ,•  l  „  ,• 


rr  -i  mi  ki  2jxl 

+  E  E  E  E[y  y.y  ]o  a  a  J) 

2=1  m=l  k=l  nrK  ] 


r  444  1 1  i  ?  k  3 

+  6  |E[w.w9w,]  e  e  e  EO-y.yJa  oJ  e> 

1  1  d  J  i=l  j=l  k=l  1  J  K 

4  4  4  .. 

+  E[w,w9w.]  e  e  e  E[y .y -y. ]o  oj  o 
1  1  4  i=l  j=l  k=l  1  J  k 

4  A  4  i  1  \  T,  k  4 

+  Ef.w-.w-w.]  e  ee  E[y.y.y. ]o  oJ  a 

1  J  4  i=l  j=l  k*l  1  J  K 


4  4  4  i  ?  •>  3  kdi 

+  E[w-w_w.]  E  E  E  E[y.y.y.]a  oJ  o  } 
L  6  4  i-1  j=l  k=l  1  J  k  ‘ 
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APPENDIX  VIII 

Generation  of  Univariate  Random  Samples 

Most  computer  systems  only  supply  random  numbers  from  the  Uniform 
distributions:  U[0,1 ).  The  best  of  these  random  number  generators  is 
discussed  in  Appendix  IX.  To  generate  data  from  multivariate  distributions, 
it  is  necessary  to  use  observations  from  other  univariate  distributions. 

This  appendix  presents  methods  to  generate  data  from  the  following  distributions: 
U(a,b),  E(A),w(b,c),  N(u,o),  LN(y.B).  and  C(a,e). 

UNIFORM 

The  uniform  distribution  has  a  distribution  function 


U(a,b) 


A  standard  technique  used  for  many  distributions  is  Inversion.  Recall  that 
any  distribution,  r  *  F(x),  is  distributed  U[0,1  ).  Therefore  if  F  is 
invertible,  x  *  F”^(r)  will  be  from  the  desired  distribution.  Specifically 

x  =  a  +  (b-a)r 


EXPONENTIAL 


The  exponential  distribution  also  uses  the  Inversion  technique  to 
generate  a  random  observation  given  an  observation  from  a  U(0,1).  The 
exponential  is 


E( A )  *  1  -  exp(-Ax) 


Set  r-  F(a).  Since  r  is  random  between  0  and  1,  so  is  ltr.  Thus 

r*  =  i-r  s  exp(-A  x) 
ln(r*)  *  -Ax 

x  =  -Mr*) 

A 


*  »*  »  •  .  • 

••>>>: 
•-  .\.V  '• 

••  y.-; 
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WE I BULL 


The  Welbull  distribution  Is  given  by 


W(b.c)  *  1  -  exp[-(§  ] 


By  the  same  analysis  as  before 


x  =  b(-  lnr)1/c 


NORMAL 


Some  distributions  cannot  be  inverted  easily.  These  must  be  handled 


by  other  methods.  The  Normal  distribution  is  the  most  Important  in  our 


study.  Let  r^,  and  r 2  be  observations  from  a  U[0,1). 


x1  *  o*c7Tnrjcos(2irr2)  +  u 


x2  *  o/^ln^sinUirrg)  +  v 


are  two  observations  from  a  N  (p,o).  (Note  that  the  operators  sin  and 
cos  act  on  radian  arguments.  ) 


LOGNORMAL 


The  Lognormal  distribution  can  be  generated  easily  once  we  have 


two  normal  observations.  Thus,  to  generate  from  a  Ln(Y»e)  ,  set 


u  *  lny  -  Jjlnf-^-  +  1) 
Y 


1  ln(  -y-  ♦  1  ) 
Y 


«.* 


1 


Then,  if  n1 .n2~N(iu o) 


are  from  the  desired  distribution. 

GAMMA 

The  gamma  distribution  Is  generated  by  summing  observations  from  the 
exponential  distribution.  Specifically 
r 

l  EUMS(X.r) 


CAUCHY 

The  Cauchy  is  a  symmetric,  fat  tailed  distribution  with  parameters  median 
a  and  scale  6.  c  gives  a  measure  of  location  while  e  provides  an  indication 
of  dispersion.  The  probability  density  function  for  the  Cauchy  is 

C(a.B)  -i  Lh(*-^-)Y’ 

The  following  generation  techniques  will  lead  to  the  generation  of  C(r,q), 
where  r  Is  real  and  q  is  rational. 

Let  Nj  and  N2  be  independently  Identically  distributed  N(0,1).  Then 


Next,  use  the  fact  that 


APPENDIX  IX 


An  Evaluation  of  Several  Random  Number  Generators 

When  we  began  performing  the  simulations  of  normally  distributed  data, 
some  anomalies  were  present  in  our  results.  Consequently,  it  was  decided 
to  test  the  quality  of  the  random  number  generator  used  in  the  simulation 
program,  and  determine  if  there  were  any  alternative  generators  which 
performed  better. 

The  following  random  number  generators  were  studied: 

1)  The  standard  VAX  BASIC  random  number  generator. 

2)  MTH-RANDOM  a  VAX  Run-Time  Library  procedure. 

3)  GGUW,  the  IMSL  routine  for  generating  random  numbers  with 
shuffling. 

4)  The  algorithm  RN0=(25173*RND+13849)  MOD  65536,  which  is  included 
in  Peter  Grogono's  book  PROGRAMMING  IN  PASCAL. 

5)  The  algorithm  RND=(1061*RND+9533)  MOD  65536,  which  was  developed 
by  G.  Silberberg  for  use  as  part  of  a  previous  project. 

A  Kolmogorov- Smirnov  test  was  used  to  measure  the  randomness  of  the 
generators.  The  Kolmogorov -Smirnov  test  is  a  standard  statistical  procedure 
to  determine  whether  a  set  of  data  can  be  generated  from  a  specified 
distribution.  The  K-S  statistic  D  is  defined  as 

D  «  max  (F(1/n)  -  FQ(1),  Fq(1)  -  F((i-l)/n)} 

where  FQ(1/n)  is  the  1th  ordered  observation  in  the  data  set  and  F  is 
the  distribution  function  which  the  data  is  to  be  tested  against. 

If  D  is  greater  the  1.22//"n  (where  n  is  the  number  of  observations 
in  the  data  set),  then  it  can  be  stated  with  90%  confidence,  that  the  data 


does  not  follow  the  given  distribution. 


The  following  results  were  obtained  using  a 

50  trials  were  performed  for  each  algorithm. 

sample  of  199  numbers 

#  Trials 

Generator 

Average  D 

Randomness  Rejected 

BASIC 

.06143 

4 

VAX  RTL 

.06308 

8 

INSL 

.06058 

5 

Grogono's 

.25837 

50 

Silberberg's 

.07324 

10 

Another  set  of  trials  was  performed,  this  time  with  a  sample  size  of 

2000. 


#  Trials  Rejected 


Generator 

Average  0 

out  of  20 

BASIC 

.01269 

0 

VAX  RTL 

.01419 

0 

IMSL 

.02314 

4 

Grogono's 

.31935 

20 

Silberberg's 

.07074 

20 

Conclusions:  The  BASIC  random  nunber  generator  showed  the  best 


performance,  and  was  placed  in  our  simulation  program.  Unfortunately,  we  found 
problems  in  sending  the  random  numbers  from  BASIC  to  Pascal.  The  system 
generator  used  by  BASIC  Is  MTHSRANDOM;  therefore,  we  used  this  system  function 
to  generate  random  numbers  with  the  Pascal  program.  The  RTL  generator  performed 
nearly  as  well,  and  may  be  appropriate  for  use  in  certain  circumstances.  Th^ 


others,  especially  Grogono's  algorithm,  should  be  avoided. 
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